function log_like = LLcorr_or(param)

%% GLOBAL

global Y T1 T2 SM SR Draws M Xc Gmax Gmin or orm ll


%% PARAM def 
cm = param(1);
cr = param(2);
phir = param(3);
c1 = param(4);
phi1 = param(5);
sigem = param(6);
siger = param(7);
sigtheta = param(8);
c2 = param(9);
phi2 = param(10);
cy = param(11);
phim = param(12);
gamma1 = param(13);
gamma2 = param(14);
by = param(15);
b1 = param(16);
b2 = param(17);
cgpa = param(18);
phigpa = param(19);
sigegpa = param(20);

c1or = param(21);
phi1or = param(22);
c2or = param(23);
phi2or = param(24);
gamma1or = param(25);
gamma2or = param(26);
% byor = param(26);
b1or = param(27);
b2or = param(28);

sige12 = param(29);

sige1 = 1;
sige2 = 1;
sigy = 1;

%%
% FunctionNo  1. logit 2. probit 3. linear

theta = Draws*sigtheta;

Lsm = normalpdf(repmat(SM,1,M)-cm-phim*theta,0,sigem);
Lsr = normalpdf(repmat(SR,1,M)-cr-phir*theta,0,siger);



if abs(sige12)>1
    sige12=sige12/(abs(sige12)+0.0001);
end;
Tt1=repmat(T1,1,M); Tt2=repmat(T2,1,M);

w=[reshape((2*Tt1-1).*(c1+c1or*or+phi1*theta+phi1or*or.*theta+b1*repmat(Xc,1,M)+b1or*or.*repmat(Xc,1,M)),[],1),reshape((2*Tt2-1).*(c2+c2or*orm+phi2*theta+phi2or*orm.*theta+b2*repmat(Xc,1,M)+b2or*orm.*repmat(Xc,1,M)),[],1)];
t1=reshape(Tt1,[],1); t2=reshape(Tt2,[],1);
Prob=((t1.*t2).*mvncdf(w,[0,0],[1,sige12;sige12,1]))...
    +((t1.*(1-t2)).*mvncdf(w,[0,0],[1,-sige12;-sige12,1]))...
    +(((1-t1).*t2).*mvncdf(w,[0,0],[1,-sige12;-sige12,1]))...
    +(((1-t1).*(1-t2)).*mvncdf(w,[0,0],[1,sige12;sige12,1]));

Lt12=reshape(Prob,[],M);
Lgpa = normcdf(cgpa+phigpa*theta-Gmax,0,sigegpa).*(repmat(Xc,1,M)>=Gmax)+(repmat(Xc,1,M)<Gmax & repmat(Xc,1,M)>Gmin).*normalpdf(repmat(Xc,1,M)-cgpa-phigpa*theta,0,sigegpa)+(repmat(Xc,1,M)==Gmin).*(1-normcdf(cgpa+phigpa*theta-Gmin,0,sigegpa));

%% Bias

yind = normcdf(cy+theta+by*repmat(Xc,1,M));
yind = yind+0;
% t1pr = normcdf(c1+c1or*or+phi1*theta+phi1or*or.*theta+b1*repmat(Xc,1,M)+b1or*or.*repmat(Xc,1,M));
% t2pr = normcdf(c2+c2or*orm+phi2*theta+phi2or*orm.*theta+b2*repmat(Xc,1,M)+b2or*orm.*repmat(Xc,1,M));
bias1 = repmat(T1,1,M)-yind;
bias2 = repmat(T2,1,M)-yind;

% b = (bias1+bias2);

% Y

Ly = (normcdf(cy+theta+gamma1*bias1+gamma1or*or.*bias1+gamma2*bias2+gamma2or*orm.*bias2+by*repmat(Xc,1,M), 0, sigy).^repmat(Y,1,M)).*((1-normcdf(cy+theta+gamma1*bias1+gamma1or*or.*bias1+gamma2*bias2+gamma2or*orm.*bias2+by*repmat(Xc,1,M), 0, sigy)).^(1-repmat(Y,1,M)));

%% Individual Likelihood Contribution
Lm=(((Lsm.*Lsr).*Lt12).*Ly).*Lgpa;
% Lm=((((Lsm.*Lsr).*Lt1).*Lt2).*Ly);


%% log likelihood
ll = log(mean(Lm,2));
log_like = -sum(log(mean(Lm,2)));